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Abstract 

Metastasis is a challenge in the management of renal cell carcinoma (RCC). 
The clinical implications of RCC metastasis are significant, with patients 
experiencing a poorer response to treatment and a reduced life expectancy. 
Currently, there is no curative treatment for RCC metastasis, and systemic 
chemotherapy and immunotherapy remain mainstay therapies. However, 
emerging evidence suggests that targeted therapies may provide promising 
treatment options for improving the prognosis and quality of life of patients 
with metastatic RCC. The current study aimed to investigate the underlying 
molecular mechanisms of RCC metastasis and identify key therapeutic targets 
using computational approaches. Also, a library of FDA approved drugs was 
screened against the target genes to obtain potential inhibitors that can be used 
in therapies. Differential gene expression analysis followed by functional 
enrichment and protein-protein interaction analyses led to the identification of 
CCND2 and MMP9 as key genes involved in multiple processes in metastatic 
tumors in the case of both clear cell and papillary cell renal carcinoma. Docking 
studies revealed good binding of drugs mysoline and ethisterone with CCND2 
and canrenone and evodiamine with MMP9. Molecular Dynamics simulations 
showed stable and strong binding of these drugs with CCND2, but the 
identified drugs failed to have stable interactions with MMP9. Hence, our study 
reveals mysoline and ethisterone to be potential inhibitors that target CCND2 
to control metastasis of RCC. 
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Introduction 


Renal cell carcinoma (RCC) accounts for the vast 
majority (approximately 90%), of kidney cancer 
cases. The two most prevalent subtypes of RCC are 
clear cell renal cell carcinoma (ccRCC) and papillary 
renal cell carcinoma (pRCC), with ccRCC accounting 
for roughly 70-80% of cases and pRCC for 15% of 
cases. The effective management and choice of 
treatments depend on the differentiation of these 
subgroups. Up to 30% of patients with RCC present 
with metastatic disease at the time of diagnosis, and 
up to 40% of those who undergo surgical resection 
develop metastases within 5 years, despite advances 
in RCC detection and therapy [1-3]. 

Invasion, migration, intravasation, survival in the 
bloodstream, extravasation, and colonization at 
distant sites are all steps in the complex process of 
metastasis [4]. The molecular mechanisms underlying 
the metastasis of RCC have been the subject of several 
investigations, all with the aim of identifying new 
treatment targets and enhancing patient outcomes. 
The hypoxia-inducible factor (HIF) pathway, the 
phosphatidylinositol 3-kinase (PI3K) pathway, and 
the mitogen-activated protein kinase (MAPK) 
pathway are three important signaling pathways that 
have been linked to the metastasis of RCC [5]. 
Moreover, many genes and proteins, such as vascular 
endothelial growth factor (VEGF), carbonic 
anhydrase IX (CAIX), and matrix metalloproteinases 
(MMPs), have been discovered as possible 
biomarkers for predicting metastatic risk in RCC ([6, 
7]). Recently developed targeted treatments have 
demonstrated promising outcomes for advanced RCC. 
The survival percentage of patients with advanced 
RCC is still dismal, despite recent advancements in 
cancer treatment. New treatment approaches are thus 
required that can successfully target the molecular 
processes underlying RCC. 

Drug research and design have made extensive use of 
computational methods like docking and molecular 
dynamics (MD) simulations. Numerous research has 
indicated encouraging outcomes when developing 
novel RCC therapeutics utilising MD simulations and 
docking. For instance, mTOR, a crucial regulator of 
cell growth and proliferation, was the subject of a 
study by Karaman et al. (2020) that employed docking 
and MD simulations to find potential inhibitors [8]. 
The investigation discovered several drugs that 
demonstrated encouraging mTOR inhibitory activity 
and may one day be used to create novel RCC 
treatments. Similarly, Sözen et al. (2019) investigated 
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the interactions of possible medication candidates 
with the VEGF receptor, a crucial regulator of 
angiogenesis in RCC, using docking and MD 
simulations [9]. The study found inhibitors with a 
strong affinity for the VEGF receptor and showed 
both in vitro and in vivo evidence of their anti- 
angiogenic activity. Docking and MD simulations 
have shown great promise in the development of new 
therapies for RCC. 

In this study, employing various genomic analysis, we 
were able to identify the key pathways and biological 
processes enriched only in metastatic RCC patients 
and identified key genes namely CCND2 and MMP9 
involved in this metastatic progression. Using virtual 
screening, we identified top drugs that dock well to the 
candidate genes CCND2 and MMP9. We further 
employed Molecular Dynamics simulations to 
validate the binding of the drugs with the candidate 
genes. We conclude that these identified potential 
drug candidates have significant binding affinity and 
selectivity for their target proteins and these can be 
further validated with experiments for the 
development of more effective therapies for RCC 
patients. 


Materials and Methods 
Data preprocessing and DEG screening 


The gene expression profile dataset GSE43477 was 
obtained from the Gene Expression Omnibus (GEO) 
database (https:/Awww.ncbi.nlm.nih.gov/geo/). This 
dataset, contributed by Narimatsu et al., comprises 20 
primary clear cell renal cell carcinomas (ccRCCs) and 
their corresponding 20 metastases [10]. Probes not 
associated with gene symbols were excluded after 
annotation. The analysis focused on identifying 
Differentially Expressed Genes (DEGs) between 
patients with metastatic clear cell renal cell carcinoma 
and those with primary clear cell renal cell carcinoma. 
A cut-off threshold was applied, requiring a -log10 P- 
value greater than 5 and a log2 fold change exceeding 
0.75 or falling below -0.75. Additionally, the gene 
expression profile dataset GSE85258, specific to 
papillary cell renal cell carcinoma, was retrieved from 
the Gene Expression Omnibus (GEO) database. 


Functional and pathway enrichment analysis 
Functional enrichment analysis is employed to 


classify functions of large-scale genomic or 
transcriptomic data according to three criteria: 


15 


Molecular function (MF), cellular component (CC), 
and biological process (BP). DAVID, which stands for 
the Database for Annotation, Visualization, and 
Integrated Discovery  (http://david.nciferf.gow/), 
functions as a comprehensive biological knowledge 
repository equipped with analytical tools designed for 
the systematic and integrated examination of 
extensive gene lists [11]. In this study, Gene Ontology 
(GO) terms and KEGG signaling pathway enrichment 
analyses for Differentially Expressed Genes (DEGs) 
were conducted using DAVID. The cut-off threshold 
was set at P<0.1. 


Protein-protein interaction network 


The Search Tool for the Retrieval of Interacting 
Genes/Proteins (STRING; http://www.string-db.org/) 
is an online database that evaluates and integrates 
Protein-Protein Interactions (PPIs). Significantly 
upregulated genes serve as queries to extract protein- 
protein interaction networks [12]. Crosstalk genes 
were identified based on filtering criteria, requiring 
association with at least three signaling pathways. 


Overall survival and expression analysis 


Overall survival analyses for key candidate DEGs 
were performed using TCGA XENA 
(http://xena.ucsc.edu/) [13]. 


Docking 


Both CCND2 and MMP9 are not crystalized for their 
entire structure. Hence, we used recently developed 
Alphafold database to obtain the models of CCND2 
and MMP9 respectively [14]. Both models underwent 
preprocessing steps, including the addition of polar 
hydrogen bonds and Kollman charges, followed by 
the generation of pdbqt structures using 
Autodocktools [15]. Subsequently, grid boxes were 
defined for the models using Autodocktools. A ligand 
library consisting of 1000 FDA-approved drugs was 
energy minimized and converted to pdbqt format 
using the PyRx tool [15]. Virtual screening was 
conducted using Autodock Vina [16, 17]. 


Molecular Dynamics Simulation 


Each of the four docked structures was utilized as the 
initial conformation for solvation employing the 
TIP3P water model within a cubic box. Standard 
minimization and equilibration methods were 
employed to bring all systems into equilibrium. 
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GROMACS [18] and the CHARMM all-atom force 
field were used for simulations. A real space cut-off 
distance of 1.2 nm was selected, and periodic 
boundary conditions were applied. For computing 
forces and potentials between grid points, a hybrid 
method was used, combining fourth-order cubic 
interpolation with particle mesh Ewald (PME) 
summation, utilizing a grid spacing of 0.16 nm. A van 
der Waals cutoff distance of 1 nanometer was applied. 
Numerical integration of the equations of motion was 
conducted using a time step of 2fs, and coordinates 
were recorded every 100 ps. The initial systems 
underwent energy minimization using the steepest 
descent approach. Prior to commencing the 100 ns 
production run, a modified Berendsen thermostat was 
employed to maintain a temperature of 310K and a 
Parrinello-Rahman barostat was used for pressure 
coupling at 1 bar. 


Analysis of Trajectories 


The global analysis of trajectories was conducted by 
calculating Root Mean Square Deviation (RMSD) 
using GROMACS modules rms. The number of 
hydrogen bonds was determined using the hbond 
module. The binding energy of the RNA with the 
protein was computed using the g mmpbsa package 
[19]. Plotly was utilized to construct graphs and plots. 
VMD was employed for visualization and the 
generation of molecular images [20]. Volcano plots 
and Dot plots were generated in RStudio [20]. using 
ggplot. High-quality images were produced using 
GIMP (GNU Image Manipulation Program) v2.10. 


Results 


DEGs identification and Functional enrichment 
analysis 


The microarray dataset GSE43477 which included 20 
primary ccRCCs and 20 corresponding metastases 
samples was taken to identify the genes and pathways 
enriched only in the metastatic tumours. Differential 
gene expression analysis was performed and the 
significantly enriched genes with log2 (foldchange) 
cut-off of 0.75 and —log10 (p-value) cut-off of 5 were 
identified (Fig. la). 78 genes were significantly 
upregulated in metastatic ccRCC. These genes were 
used for functional enrichment using DAVID. The 
KEGG pathways, biological processes and cellular 
components enriched by these genes are identified 
(Fig. 1b-d). TGF-beta signalling pathway, 
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Fig. 1: Differential expression and functional enrichment of genes upregulated in metastatic patients reveals major 
signalling pathways in the disease. (a) Volcano plot depicting the differentially expressed genes between metastatic and primary 


clear cell renal carcinoma patients. Cut-offs: log2(foldchange) = 


0.75 and -log10 (p-value) = 5. (b) Dot plot depicting the KEGG 


pathways enriched by significantly upregulated genes in metastatic clear cell renal carcinoma patients. P-value cutoff = 0.05. (c) 
Dot plot depicting the biological processes enriched by significantly upregulated genes in metastatic clear cell renal carcinoma 
patients. P-value cutoff = 0.05. (d) Dot plot depicting the cellular components enriched by significantly upregulated genes in 
metastatic clear cell renal carcinoma patients. P-value cutoff = 0.05 


parathyroid hormone synthesis, oocyte meiosis, 
human T-cell leukaemia virus infection, hippo 
signalling pathway and estrogen signalling pathway 
are the major pathways enriched by the upregulated 
genes. BMP signalling process, apoptotic process and 
regulation of inflammatory response are the processes 
significantly enriched. Cell membrane, cytosol and 
plasma membrane are the cellular components 
enriched by the genes upregulated in the metastatic 
ccRCC. 


Crosstalk genes and key candidate gene screening 


The significant upregulated 78 genes are then 
subjected to string network analysis to identify the 
protein-protein interaction networks among them 
(Fig. 2a). The genes that are enriched in at least 3 
different processes and have at least 2 interactions 
with other proteins are taken as the candidate subjects. 
BMP7, CCND2, EGR1, MAPK3, MMP9 are the 
candidate genes identified through this process. The 
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survival rate of the ccRCC patients with differential 
expression of these candidate genes are analysed using 
the TCGA dataset. CCND2 and MMP9 are the only 
genes that resulted in lower survival probability of the 
ccRCC patients when their expression is high (Fig. 
2b-f). 

We also checked the expression of CCND2 and 
MMP9 in papillary renal cell carcinomas using 
GSE85258 microarray dataset and identified their 
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expression to be high in metastatic carcinomas 
compared to primary carcinomas (Fig 3a). The 
survival probability of pRCC patients is also 
significantly lower when the expression of CCND2 
and MMP9 is high. From these analyses, we conclude 
that CCND2 and MMP9 are potential candidates that 
can be used as a therapeutic target to stop the 
progression of RCC to metastasis. 
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Fig. 2: Network and survival analysis reveals the important genes for therapeutics. (a) Network analysis using string database 
for the significantly upregulated genes. (b-f) Survival probability of clear cell renal carcinoma patients in TCGA dataset with 
differential expression in the genes BMP7, CCND2, EGR1, MAPK3, and MMP9, respectively. 
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Fig. 3: Gene Expression and survival analysis in papillary renal cell carcinoma patients in TCGA dataset also reveals MMP9 
and CCN2 as candidate genes. (a) Gene expression of MMP9 and CCND2 between Metastatic and Primary papillary renal cell 
carcinoma patients. (b) Survival probability of papillary cell renal carcinoma patients in TCGA dataset with differential expression 


in the genes CCND2 and MMP% respectively. 


Virtual Screening to 
approved drugs 


identify potential FDA 


Molecular docking is a useful method for analysing a 
compound's potential for binding to a certain target 
protein. Many now employ this technique in their 
efforts to find new drugs, which has gained 
widespread recognition. This docking output can be 
used to separate target-specific ligands with high 
potency from those with poor binding affinity. Low 
binding affinity is correlated with lower binding 
energy, and high binding affinity is correlated with 
higher binding energy. 

A library of 1000 FDA approved drugs are prepared 
and docking of this drug library with CCND2 and 


MMP9 were performed using Autodock vina. The top 
5 drugs with high binding affinity are identified. The 
number of hydrogen bonds formed by these drugs 
with CCND2 and MMP9 were also identified (Fig. 
4a). The drugs mysoline and ethisterone are the drugs 
that had a high binding affinity (-8.3 Kcal/mol and — 
7.9 Keal/mol respectively) and formed hydrogen 
bonds with CCND2 protein. Canrenone and 
evodiamine are the drugs that had a high binding 
affinity (-9.4 Kcal/mol and -8.9 Kcal/mol 
respectively) and formed hydrogen bonds with MMP9 
protein. The best-docked conformer of mysoline and 
ethisterone with CCND2 (Fig. 4b, 4c) and the best- 
docked conformer of canrenone and evodiamine with 
MMP9 (Fig 4d, 4e) were taken and used as starting 
structure for Molecular Dynamics simulations. 
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Fig. 4: Virtual screening identifies candidate FDA approved drugs to inhibit CCND2 and MMP9. a) Table representing the 
top five drugs identified by virtual screening. The binding energy of each drug with CCND2 and MMP9 and the number of 
hydrogen bonds they form are represented. b), c), d) & e) Snapshots depicting the docking of mysoline and ethisterone with 


CCND2 and canrenone and evodiamine with MMP9 respectively. 


Molecular Dynamics simulation validates the 
observations identified through virtual screening 


MD simulation studies are employed to evaluate the 
conformational stability of protein-ligand complexes. 
The best-docked conformer of mysoline and 
ethisterone with CCND2 was used as a starting 
structure and Molecular Dynamics simulation was 
performed for 100ns. Snapshots at different 
timepoints are depicted (Fig. 5a, 5b). Analysis of MD 
simulation snapshots indicates that throughout the 100 
ns simulations, each ligand interacts with the target 
protein through various binding poses. Root mean 


square deviations (RMSD) were calculated for the 
protein CCND2 to assess its stability in the presence 
of the two ligands. The RMSD plot demonstrates that 
protein CCND2 remained stable in the presence of 
both mysoline and ethisterone throughout the 
simulation, with average RMSD values of 0.479 nm 
and 0.478 nm, respectively (Fig. 5c). The stability of 
the initially docked ligand at the protein's active 
binding site is determined by the ligand RMSD. 
Ligand RMSD little higher than the protein RMSD is 
permissible. When a ligand's RMSD varies a lot, it 
indicates that it is striving to find a new pose as the 
simulation progresses. Ligand RMSD for both 
mysoline and ethisterone are stable throughout the 
simulation with an average RMSD of 0.024 nm and 
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0.14 nm respectively (Fig. 5d). MM-PBSA analysis 
was employed to estimate the binding free energies of 
protein-ligand complexes, aiming for low values, 
which indicate potent drug candidates. MM-PBSA 
offers greater precision and reliability compared to 
molecular docking scoring functions for calculating 
binding energies. To enhance the accuracy of 
assessing the binding strength and potency of 
mysoline and ethisterone against CCND2, binding 
free energies for all protein-ligand complexes were 
computed throughout the 100 ns MD simulation. 
Average binding energy of mysoline and ethisterone 
with CCND2 are —64.17 KJ/mol and —-42.34 KJ/mol 
respectively (Fig. 5e). During the simulation, it's 
crucial to observe protein-ligand interactions to gauge 
the ligand's capability to interact with amino acids in 
the protein. Therefore, hydrogen bond interactions 
between the ligand and protein are computed 
continuously throughout the simulations to gain 
insights into their interactions (Fig. 5f). We observe 
that mysoline had more hydrogen bonds formations 
with CCND2 compared to ethisterone. From these 
analyses, we conclude mysoline to be better potent 
inhibitor of CCND2 as it had better binding potential 
and more hydrogen bond interactions and can be taken 
ahead for experimental validations. 

The best docked conformer of canrenone and 
evodiamine with MMP9 was used as a starting 
structure and Molecular Dynamics simulation was 
performed for 100ns. Snapshots at different 
timepoints are depicted (Fig. 6a, 6b). Analysis of MD 
simulation snapshots reveals that over the entire 100 
ns duration, each ligand adopts various binding poses 
while interacting with the target protein. Root mean 
square deviations (RMSD) were computed for the 
protein MMP9 to assess its stability in the presence of 
the two ligands. RMSD plot revealed that the protein 
MMP9 was stable after a few initial fluctuations in the 
ligand bound state throughout the simulation with an 
average RMSD of 0.91 nm and 1.1 nm respectively 
for canrenone and evodiamine bound simulations 
(Fig. 6c). The protein underwent approximately Inm 
deviation in the first 10ns and then remained stable. 
Ligand RMSD for both canrenone and evodiamine are 
calculated. Canrenone was stable throughout the 
simulation with an average RMSD of 0.03 nm while 
evodiamine fluctuated a lot stable throughout the 
simulation (Fig. 6d). To get a more accurate 
assessment of the binding strength and potency of 
canrenone and evodiamine against MMP9, the 
binding free energies for all protein-ligand complexes 
were calculated during the 100 ns MD simulation. 
Average binding energy of canrenone and evodiamine 
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with MMP9 are —75.57 KJ/mol and —72.23 KJ/mol 
respectively (Fig. 6e). The hydrogen bond interactions 
formed between the ligand and protein are calculated 
throughout the simulations (Fig. 6f). We observe that 
evodiamine lost its hydrogen bonds with MMP9 in the 
beginning of the simulation and did not form any 
hydrogen bonds. Whereas canrenone had formed 1 
hydrogen bond at 5 different timepoints but a stable 
interaction was not observed. From these analyses, 
we conclude that even though both the ligands have 
good binding affinity, they lack any hydrogen bond 
interaction with the MMP9 protein. The best potent 
drugs identified using docking, failed to have stable 
interactions in dynamics conditions and thus do not 
qualify for further experimental validations. 


Discussion 


The intricacy of the RCC disease's management and 
therapy become much more difficult when RCC 
metastasize. Determining possible therapeutic targets 
requires an understanding of the mechanisms behind 
metastasis. A number of intricate processes are 
involved in the formation of metastases, including 
Local Invasion, Intravasation, Circulation, Arrest and 
Extravasation, Colonization and Formation of 
Metastases. Several strategies were taken into 
consideration to successfully target metastasis in 
RCC. The development of new blood vessels, or 
angiogenesis, is essential for the growth and 
metastasis of tumors. Angiogenesis-inhibiting drugs, 
such as sunitinib and bevacizumab, have been 
employed in the treatment of metastatic RCC to 
disrupt the tumor's blood supply and hinder its 
dissemination [21, 22], Similarly, tyrosine kinase 
inhibitors (TKIs) target several enzymes involved in 
cell signaling pathways that promote tumor growth 
and metastasis. Given the frequent overactivity of 
tyrosine kinases in RCC, medications like sorafenib, 
pazopanib, and axitinib work to reduce their activity, 
thereby impeding the progression and spread of the 
malignancy [23]. Immunotherapy strengthens the 
immune system's defences against cancerous cells. 
Immune checkpoint inhibitors, like ipilimumab and 
nivolumab, can assist in directing the immune 
system's attack toward and eliminating cancer cells, 
even those that have spread [24]. The identification of 
certain molecules and pathways involved in the RCC- 
specific metastatic process is still under investigation. 
It could be able to stop or reduce metastasis by 
focusing on these pathways. In clinical trials, HIF-20 
pathway inhibitors like belzutifan have demonstrated 
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potential in treating metastatic RCC [25]. Treatment 
outcomes can be enhanced by combining diverse 
modalities, such as immunotherapy or chemotherapy 
with targeted therapy, which works synergistically to 
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Although these treatments have produced 


encouraging outcomes, there is a need for other 
therapeutic choices because treatment resistance is 
widespread. 


target multiple features of RCC, including metastasis. 
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Fig. 5: Molecular Dynamics simulation of CCND2 with Mysoline and Ethisterone. (a) Snapshots at different time points during 
simulation of mysoline with CCND2. (b) Snapshots at different time points during simulation of ethisterone with CCND2. (c) The 
root mean square deviation of the protein CCND2 during the course of simulation is represented as line plot. Simulation of CCN2 
with mysoline and ethisterone are represented in orange and green respectively. MYS — mysoline bound CCN2 simulation, ETS — 
ethisterone bound CCND2 simulation. (d) The root mean square deviation of the drugs during the course of simulation is 
represented as line plot. Simulation of CCN2 with mysoline and ethisterone are represented in orange and green respectively. (e) 
The binding energy between the protein CCND2 and the drugs during the course of simulation is represented as line plot. Simulation 
of CCN2 with mysoline and ethisterone are represented in orange and green respectively. (f) The number of hydrogen bonds 
between the protein CCND2 and the drugs during the course of simulation is represented as bar plot. Simulation of CCN2 with 
mysoline and ethisterone are represented in orange and green, respectively. 
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Figure 6. Molecular Dynamics simulation of MMP9 with Canrenone and Evodiamine. (a) Snapshots at different time points 
during simulation of Canrenone with MMP9. (b) Snapshots at different time points during simulation of Evodiamine with MMP9. 
(c) The root mean square deviation of the protein MMP9 during the course of simulation is represented as line plot. Simulation of 
MMP9 with Canrenone and Evodiamine are represented in violet and blue respectively. CNN — Canrenone bound MMP9 
simulation, EDM- Evodiamine bound MMP9 simulation. (d) The root mean square deviation of the drugs during the course of 
simulation is represented as line plot. Simulation of MMP9 with Canrenone and Evodiamine are represented in violet and blue 
respectively. (e) The binding energy between the protein CCND2 and the drugs during the course of simulation is represented as 
line plot. Simulation of MMP9 with Canrenone and Evodiamine are represented in violet and blue respectively. (£) The number of 
hydrogen bonds between the protein CCND2 and the drugs during the course of simulation is represented as bar plot. Simulation 
of MMP9 with Canrenone and Evodiamine are represented in violet and blue respectively. 


In preclinical research and clinical trials, the 
development of drugs that target important signalling 
pathways such the PI3K/Akt/mTOR, TGF-, and HIF- 
2 signalling pathways has shown promising outcomes. 


Future studies in the hunt for new therapeutic targets 
for RCC metastatic treatments should concentrate on 
creating safer, more efficient treatments. Integrating 
genomic analysis and drug discovery to identify and 
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target key genes involved in the major pathways of 
metastatic tumours is the focus of the study. With this, 
we were able to identify key genes namely CCND2 
and MMP9 enriched significantly only in metastatic 
tumour samples. Virtual screening helped us to 
identify top 2 drugs that bound with these proteins. 
MD simulations revealed drugs bound to MMP9 
failed to have stable interactions whereas drugs bound 
to CCND2 formed stable interactions throughout the 
simulation. We propose that mysoline and ethisterone 
drugs can be repurposed to treat patients with 
metastatic renal carcinomas with further experimental 
validations. 
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